%% boundary_edges (for Neumann and Robin)
function bde = boundary_edges(mesh, index)
gap_x = mesh.Nx;
gap_y = mesh.Ny;
switch index
    case 1 % bottom
        bde = mesh.bde(1:gap_x,1:3);
    case 2 % right
        bde = mesh.bde(gap_x+1:gap_x+gap_y,1:3);
    case 3 % top
        bde = mesh.bde(gap_x+gap_y+1:2*gap_x+gap_y,1:3);
    case 4 % left
        bde = mesh.bde(2*gap_x+gap_y+1:end,1:3);
    otherwise
        error("Invalid boundary number.");
end
end